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A MULTI-DOMAIN SPECTRAL METHOD FOR SUPERSONIC REACTIVE FLOWS * 


WAI-SUN DON, DAVID GOTTLIEB & JAE-HUN JUNG 

Abstract. This paper has a dual purpose: it presents a multidomain Chebyshev method for the solu- 
tion of the two-dimensional reactive compressible Navier-Stokes equations, and it reports the results of the 
application of this code to the numerical simulations of high Mach number reactive flows in recessed cavity. 
The computational method utilizes newly derived interface boundary conditions as well as an adaptive fil- 
tering technique to stabilize the computations. The results of the simulations are relevant to recessed cavity 
flameholders. 

Key words, multi-domain spectral method, penalty interface conditions, supersonic combustor, recessed 
cavity flame-holder, compressible Navier-Stokes equations 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. The efficacy of spectral methods for the numerical solution of highly supersonic, 
reactive flows had been previously reported in the literature. Don and Gottlieb [7, 8] simulated interactions 
of shock waves with hydrogen jets and obtained results showing the rich dynamics of the mixing process as 
well as the very complex shock structures. Don and Quillen [9] studied the interaction of a planar shock 
with a cylindrical volume of a light gas and showed that the spectral methods used gave good results for the 
flows with the shocks and complicated non- linear behaviors, in fact the results compared favorably to ENO 
schemes. 

The methods reported above were based on Chebyshev techniques in one domain. In order to extend 
the utility of spectral methods to complex domains, multidomain techniques have to be considered. The 
main issue here is the stable imposition of the interface boundary conditions, and in this paper we consider 
mainly the penalty method, introduced for hyperbolic equations by Funaro and Gottlieb [10, 11]. 

There is an extensive literature on the subject: Hesthaven [13, 14, 15] applied penalty BC for Chebyshev 
multidomain methods using the characteristic variables. Carpenter et. al. [4, 17, 18] used it in conjunction 
with compact finite difference schemes, going from a scalar model equation to the full N-S equations in 
general coordinate systems. Carpenter, Gottlieb and Shu [5] demonstrated the conservation properties of 
the Legendre multidomain techniques. 

Iri the current work we follow the same methodology but in the context of supersonic combustion. 
We formulate the stable interface conditions based on the penalty method in a conservative form for both 
Euler and Navier-Stokes equations in two dimensional Cartesian coordinates. We derive stability conditions, 
independent on the local flow properties, for the penalty parameters for the Legendre spectral method. 
We also present here a new adaptive filtering technique that stabilize the spectral scheme when applied to 
supersonic reactive flows. 

Implementing this method, we consider supersonic combustion problems in recessed cavities in order to 
establish the efficacy of recessed cavity flame-holders. 

* Brown University, Division of Applied Mathematics, 182 George Street, Providence, RI 02912 (E-mail: wsdon, dig, 
jung@cfm.brown.edu) This work was performed under AFOSR grant no. F49620-02-1-0113 and DOE grant no. DE-FG02- 
96ER25346. This research was supported by the National Aeronautics and Space Administration while the second author was 
in residence at ICASE, NASA Langley Research Center, Hampton, VA 23681. 
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We consider two different cases; (1) Non-reactive flown with two chemical species and (2 ) Reactive flows 
with four chemical species. 

Recessed cavities provide a high temperature, low speed recirculating region that can support the pro- 
duction of radicals created during chemical reactions. This stable and efficient flame-holding performance 
by the cavity is achieved by generating a recirculation region inside the cavity where a hot pool of radicals 
forms resulting in reducing the induction time and thus obtaining the auto-ignition [2, 22]. Experiments have 
shown that such efficiency depends on the geometry of the cavity such as the degree of the slantness of the 
aft wall and the length to depth ratio of cavity L/D. Thus one can optimize the flame-holding performance 
by properly adjusting the geometrical parameters of the cavity flame-holder system for a given supersonic 
flight regime. There are two major issues of such cavity flame-holder system that need to be investigated ; 
{l)What is the optimal angle of the aft wall for a given LjD ? and (2) How does the fuel injection interact 
with cavity flows? An answer to these questions require both a comprehensive laboratory and numerical 
experiments. 

There have been previous numerical studies on these questions, many of them rely on the turbulence 
models. Rizzetta [19] used a modification of the Baldwin Lomax algebraic turbulence model. Davis and 
Bowersox [6] also used Baldwin- Lomax model. Zhang et.al. [23] used Wilcox k-oj turbulence model. Baurle 
and Gruber [3] used the Menter model. Although the use of the turbulence models can make it possible 
to handle the compressible supersonic shear flows, the results are quite model-dependent as they require 
parametric assumptions. In this work, we solve the full compressible Navier-Stokes equations with chemical 
reactions without any turbulence model, using a multi-domain spectral method. 

Results of several numerical studies including the present study have shown that the stability of the 
recirculation inside cavity is enhanced for the lower angle of cavity compared to the rectangular cavity. 
The present study, however, gives more accurate and finer details of the fields than those done by lower 
order numerical experiments. We show that a stationary recirculation region is not formed inside the cavity 
contrary to what the lower order schemes predict. A quantitative analysis made in this study shows that 
the lower angled wall of the cavity reduces the pressure fluctuations significantly inside the cavity for the 
non-reactive flows. We obtained a similar result for the reactive flows with the ignition of the fuel supplied 
initially in the cavity. 

The rest of this paper is organized as follows. In section 2 the governing equations are given. In section 
3 we describe the numerical method used in this work. In this section we present the adaptive-filtering 
used to remove the high frequency mode that causes the instability due to the non-smoothness of the flow, 
and we derive stable penalty interface conditions. In section 4 the system of the supersonic recessed cavity 
combustor is described. In section 5 the main results of this work are given and discussed. 

2. The Governing Equations. In this work, we consider the compressible Navier-Stokes equations 
in the presence of the chemical reactions. Since Hydrogen is used as a fuel in our numerical experiments, 
four chemical species are considered, i.e. H 2 ,0 2 ,H 2 0 and N 2 with the chemical reaction between Hydrogen 
and Oxygen gases: 

(2.1) 2H 2 + 0 2 # 2H 2 0. 


The two-dimensional compressible Navier-Stokes equations in conservative form can be w r ritten as: 


( 2 . 2 ) 


dq 81- d G _dF^ 3G V 

dt + dx + 8y dx + dy 
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The state vector, q and the inviscid fluxes, F and G are given by 
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Here p is the density, u and v are the mean mixture velocity components of flow, E is the total internal 
energy and P is the pressure. The mass fraction vector, is f = (/i , f- 2 -, h- fi) 1 and the column vectors f u 
and f v are composed of the specific momentum of i th species 


(2-4) 


fu i = fi (« + Ui ) , f V i = fi (v + Vi ) . 


The velocity field (ui,Z\) of the i th species is the drift velocity relative to the mean mixture velocity (u.v) 
and is determined by 


(2.5) 


(«»,«*) 


PS, 


■v/ 4 . 


Here p is the mixture dynamic viscosity to be determined in (2.11), and S c is the Schmidt number which is 
taken to be 0.22, The viscous fluxes, F„ and G v are given by 


( 2 . 6 ) 
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where 0 = (0, 0,0, 0) T , T is the temperature, C p is the mixture specific heat at constant pressure, P T is 
the Prandtl number (which is taken to be 0.72) for the normal air and hi is the specific enthalpy of the i th 
species and given by 


hi 


= h° 4 + 



ds. 


where ft® is the reference enthalpy of the i th species and the specific heat, of the i th species at constant 
pressure, C Pi is represented as a fourth-order polynomial of T (see [16j). The elements of the viscous stress 
tensor are given by 


(2.7) 


TxiXj ft 



duj 

Ox., 



where 5 is the Kronecker delta symbol, and A is the bulk viscosity which is taken to be — under the Stokes 
hypothesis. 
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The equation of state is given by the assumption of the perfect gas law 

4 

(2.8) P = pRT = RT V pfi/Mi, 

i= 1 

where R is a mixture gas constant with the universal gas constant R and M-, is the molecular weights of i th 
species. The energy E is given by 

.rj 1 4 

(2.9) E = p f C p (s)ds -P + \p(u 2 + v 2 ) + V p/j/t-, 

z i= 1 

where the mixture specific heat at constant pressure is given by 

4 

(2.10) C p V C Pi fi/Mi. 

i= 1 

2.1. The chemical models. We use the same models as in [7], Each chemical species has different 
dynamical viscosity p, based on Sutherland’s law and we obtain the mixture viscosity p by Wilke’s law [21], 


Pi_ 

I'Oi 


T 0i 


3/2 


T 0i + Si \ 
T + S; ) ’ 


( 2 . 11 ) 


/-E 


Pifi/Mi 


Tit fj/MjPij 
(.i+[(p t /p j )if j /f i )}HM i /Ai J )iy 


Here po it T 0i and Si are constants. A modified Arrhenius Law gives the equilibrium reaction rate k e , the 
forward reaction rate kf and the backward reaction rate /►> as 


fee = A e T exp(4.60517(i?e/T - 2.915)) 
kf = A f exp(- E f /(RT)) 
kb — kf j k e . 


where the activation energy E e = 12925, Ef = 7200 and the frequency factor A e = 83.006156, Af = 5.541 x 
10 14 . 

The species are ordered as follows : (H 2 ,02,H 2 0,N2), and the law of mass action is used to find the 
net rate of change in concentration of i th species C, by the single reaction (2.1), i.e. 

Ci = 2(fe/[H 2 ] 2 [0 2 ] - fei,[H 2 0] 2 ) 

Ci = — (fe/[H 2 ] 2 [0 2 ] - kt [H 2 0] 2 ) 

C 3 = 2(kf [H 2 ] 2 [0 2 ] - fe 6 [H 2 0] 2 ) 


where [•] denoted the net rate of change in concentration. 
Finally, the chemical source term C in (2.2) is given by 


( 2 . 12 ) 


C = 


(°* 


0, 0, 0, C] Ml , C 2 M 2 , C 3 M 3 , Ci A 


M,) 


where C» is the net rate of change in concentration of i th species by the reaction. 

In Appendix C, a table of all the necessary coefficients and constants used for the reactive Navier-Stokes 
equations with species (H 2 . 0 2 , H 2 0, N 2 ) are given. 


4 



3. The Multi-domain Spectral Method, In this section we describe the two crucial components of 
the Chebyshev multi domain code used in our work, i.e. the adaptive filtering and the penalty method for 
the stable interface conditions. 

3.1. The adaptive filtering technique. It is well known that spectral methods may exhibit insta- 
bilities when applied to nonlinear equations. To stabilize the spectral scheme in an efficient way we use here 
filters to attenuate the high frequency modes of the function q N (x,t) smoothly to zero. Thus the filtered 
version of a polynomial q N is given by: 

N 

(3.1) &v(M) = a(k/N)a k (t)(p k {x), 

k=—N 

where a* is the transform coefficient and fa is the basis polynomial of order k (generally the Fourier and 
Chebyshev polynomials for a periodic and non-periodic function respectively). 

Following Vandeven [20] we define a filter function a(u) of order p > 1 as a C'°°[— 1, 1] function satisfying 

(3 2) *(0) = 1, <r(±l) — 0 , 

(0) = 0 , a <j) (±1 .) = 0 , j < p 

where orW) denotes the j - th derivative. 

It can be shown that the filtered sum (3.1) approximates the original function very well away from the 
discontinuities. A good example of filter function is the exponential filter. It is defined as 

(3.3) a(ou') — exp (— a|w| 7 ) , 

where —1 < u = k/N < 1, a = — Inc, c is the machine zero and 7 is the order of the filter. 

The exponential filter offers the flexibility of changing the order of the filter simply by specifying a 
different 7. One does not have to write a different filter for different order. Thus varying 7 with N yields 
exponential accuracy according to [20]. in the present study the sixth order global smoothing (7/ = 6) is 
used. If the order of the filter 7 is taken to be too small, say 7 < 4. the method becomes too dissipative. 

In the current application, the interaction of the aft cavity wall and the strong vortex generated by the 
shear layer flow over the cavity, creates large pressure variations near the corner of the aft cavity wall. The 
local sharp gradient can cause numerical instability arid a heavier filter is needed to prevent the development 
of oscillations in this region. This heavy filtering can be used globally and maintain the stability of the 
scheme, however this dissipates out all fine scale structures, which is highly undesirable when the resolution 
of fine scale structures is essential for the understanding of the recessed cavity flameholder systems. 

Since this is a local phenomenon, it is enough to apply a heavy filter only in points in this region. 
This Loco, l Adaptive Filtering keeps the scheme stable, without dissipating fine scale features away from this 
region. The local adaptive filtering is carried out where conditions such as q ' < q < q u are violated. Here q 
can be the mass fraction of each species /» and/or temperature T and q l and q u denote the lower and upper 
tolerance limits of q. In this work, a filtering of the order 7 = 2 or 7 = 3 is used to reduce the magnitude of 
the oscillations at those points. 

The results of this work indicate that the local adaptive filtering is applied only in a few number (in the 
range of 1 to 7) of grid points around the corner of the aft wall once in a while. 

3.2. Stable interface conditions. In this paper we use mainly the penalty type interface conditions, 
i.e. the boundary conditions are imposed only in a weak form [10, 11]. Successful penalty interface conditions 
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were constructed based on the characteristics for the Navier-Stokes equations in [13, 14, 15] and lor spectral 
method and for high-order finite difference methods in [4, 17, 18], and a conservative form of penalty interface 
conditions was proposed [5] for the Legendre spectral method. Following the same idea as those works, we 
consider two interface conditions, i.e. 

1. The averaging method, in which the interface conditions are obtained by averaging the state vectors 
of the two adjacent domains, and 

2. The Penalty method in conservative form in which the interface conditions are satisfied only in a 
weak form, leaving the approximations not necessarily continuous at the interfaces. 

In the following sections we will give the penalty interface conditions for the Euler and Navier-Stokes 
equations and also show that the averaging method is a subset of the penalty method. 


3.2,1, Conservative penalty interface conditions. Consider equation 2.2 with the inviscid part 
only, in the ^-direction in the interval —2 < x < 2, i.e., 


(3.4) 


dq cTF 
dt + dx 


= 0 . 


For simplicity, assume that we have two domains in this interval with the interface at x = 0, q* N {x, t) 
denotes the numerical solution in the left domain x < 0 and qjj(x,t) in the right domain x > 0. Note that 
the numerical solution is composed of two polynomials of different orders. The Legendre spectral penalty 
method is given by 


(3.5) 


~m + dx -s(«*(- 2 ’0) + 

TiQ N (x)[f + {q r N {0,t)) - / + (« m ( 0 ,*))] + 


t- 2 Qn{x)[J (q!w(0,t))-f (</m(0,0)], 


Ml + dI i ! ! F{ TM) 

dt dx 


- B(«#(2,t)) + 

T 3 QM(x)[f + (q{t(0,t)) - / + (« i y(0,t))] + 


t a Qm(x)[} («m(0,0)-/ (<7w(M))] 


where B is a boundary operator at the end points, i.e., x — ±2 and I : N and /{/ are the Legendre interpolation 
operators for the left and right domains respectively. . The positive and negative fluxes / + and f~ are defined 
by 

(3.6) f± = J SA ± S~ 1 dq, 

with 

fiF 

(3.7) ,4 = — = SAS- 1 . 

dq 

The Jacobian matrix A is assumed to be symmetric. A + and A - are the diagonal matrices composed of 
positive and negative eigenvalues of A respectively. Qn(x) and Qm(x) are polynomials of orders N and M 
respectively such that they are zero at all the collocation points except the interface points x = 0 (for example 
Qn(x) — , o < x < 2 where 2’ ; y (x) is the Chebyshev polynomial of degree N). The penalty 

parameters n,r 2 ,r 3 and xj are all constants. Since we are interested only in the interface conditions, we 
ignore the boundary operator B at x = ±2. Define the discrete scalar product ( p,q)N = P T (&)<7(& )<■>>*• 

Ui is the weight in the Gauss-Lobatto-Legendre quadrature formula. With the discrete product, the energy 



Eft) is defined by E(t) = (qlf{x,t),q I N (x,t))N + (x, t), (x, ()) m- The stability conditions of penalty 

parameters are given by the following theorem [5]: 

Theorem 3.1. The energy is bounded by the initial energy of the system if the following conditions are 
satisfied ; 

V T\ £ 1, 2ujj\rT-2 > 1. 2ujfij T3 < 1, 2ul [ < T \ > 1, 

(3.8) '-OpfTl — OJmTs = 1, U)ffT2 — 0Jj{fT4 = 1. 


3.2.2. The penalty method for the Euler Equations. The penalty method in the case of the 2-D 
Euler equation is given by 


(3.9) 


dq N d I N F{ q N ) 
dt dx 


dIj\G(qx) 

dy 


= n,3 Q(x,y)[f + {q N ) - f + {q.M -)] + 


T2aQ{xaj)[J (q N )-f (q.M-)}, 


where qM- is the state vector of the adjacent domain at the interface of degree M, 1 - 1 , 3 ( 75 , 4 ) denotes 71 ( 72 ) 
and 73 ( 7 - 4 ) respectively. r\ and 75 ( 7-3 and 74 ) are the penalty parameters for the right(left) in ^-direction 
and top(bottom) in y-direction respectively. Qix.y) is a polynomial which vanishes at all of interior points 
of the domain and is equal to 1 at the four interfaces. Note that the boundary operator B does not appear 
in the scheme. Let .4 be the linearized Jacobian matrix (around a state vector q 0 ) of two inviscid fluxes 


(3.10) 


4 = 




where ft = (n x ,n v ) is the unit outward normal vector. Since the matrix A is symmetric, there exists S such 
that 


(3.11) A - S:\S~ 1 , 

where A is a diagonal matrix composed of eigenvalues of 4. Then 4 = 4 + + A~ and A ± — SA ± S~ 1 . A* is 
defined as in previous section. Splitting 4 yields 

(3.12) f ± = A ± q 0 ,. 


where f ± is obtained from the linearized state. 

Remark 1. Since it = ( n x ,n y ) is taken to be outward normal vector, the stability condition (3.8) is now 
modified and given as 


(3.13) 


2oj(yTi < 1 , 201^72 > 1 , < 1 , 2 ( 0^74 > 1 , 

OJftTi + iu'^fT4 = 1 , Lo(\jT-2 + ~ 1 . 


The Jacobian matrix A and its eigenvalue matrix A are given in Appendix A. 

For illustration, we consider the propagation of a Gaussian density peak at the center of rectangular 
physical domains. The physical domain is partitioned with 16 sub-domains. The interface conditions be- 
tween the domains are imposed according to the penalty Euler equations as discussed above. Characteristic 
boundary conditions are imposed at the outer physical boundaries. The results presented in figure 3.1 indi- 
cate that the penalty formulation works well. From the numerical experiments of this problem, we observe 





Fig. 3.1. The propagation of a density peak with the penalty Euler equations with 16 sub-domains: The initial condition 
(left) and the solution (right) at t — 0.0 3604 ms are given. 


that reflections can be created at the interface across the adjacent, domains depending on the choice of the 
penalty parameters. Thus proper choice of the penalty parameters should take into account reflections from 
the interfaces. To demonstrate the above formulation for the Euler Equations, We will return to this issue 
in a future paper 1 . 


3.2.3. The penalty method for the Navier-Stokes Equations. When dealing with the Navier 
Stokes equation, we keep the penalty form for the Euler fluxes and add a penalty term for the viscous 
fluxes. The stability of this procedure stems from the fact that the Jacobian matrices for the full reactive 
Navier-Stokes equation can be symmetrized by the same similarity transformation (see Appendix B). Thus 
we get the system: 

dq N dI N F dI N G _ cU N F„ d! N G„ 

dt * dx * dy dx + dy 

n,3 Q(x,y)[f + (qN) - / + (<zm -) ] + 
r- 2 AQ{x,y)[f~{q ,v) - f~(q M -)] + 
re. s Q(x,y)[A v ■ q A - A„ • q<w-] + 

(3.14) t 5J Q(x, y)[ A v - <9q, v - A„ ■ dq M ._]. 


Here are same as defined in the previous section and the Jacobian matrix vector A„ is given by 


( 3 . 15 ) 

and 


dF v dG„ 




dq x ' dq. 


( 3 . 16 ) 


q =(q,q), dq=(q x ,q v ), 


where again q_ and dq~ denote the adjacent domains state vectors and their derivatives. Note that the 
penalty terms A„ • dq does not appear in [4, 17, 18]. The penalty parameters tsj and re. 8 are defined in 
the same way as in the previous section. To seek stable penalty parameters we split the in viscid and viscous 
fluxes and keep the stability conditions of ti. 2 , 3.4 for the inviscid flux as in Theorem 3.1. The stability 
conditions of rg,r and rg.s are given in the following Theorem : 

1 J. H. Jung, PhD thesis, Div. of Applied Math., Brown University, 2002 
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Theorem 3.2. The penalty method for the Navier- Stokes equations (3.14) is stable if the penalty pa- 
rameters Tj , j = 1...4 are as in Theorem 3.1 and the rest satisfy: 

O \ 7(i o' 0. 

co N r e — ojmt~s = 0 , 

1 + Ou'nTs — COM r 7 = 0 , 

1 \ 9 , „ , 1 

H ) — 2t 7 + 4wj\rT6 + < 0 . 

CON J ‘ COM 


(3.11 


1 

COM 


Proof. 

As in the proof of Theorem 3.1, we assume that we have two domains and by multiplying the equations 
by the state vectors, we get 

(3.18) < [Inviscid] + [' Viscous ], 


where [Inviscid\ and [Viscous] denote the terms from inviscid and viscous parts of the equation respectively. 
The conditions for t \. 2 and 7 - 3,4 given in the Theorem 3.1 assure that the first term [Inviscid] is negative. 
The [ Viscous ] part at the interface is given by 


N 

[Viscous] = q T A„q' - ]T q [ 1 A^to; - 
i = 0 
M 

q-A„q'~ - X'/'.. 7 -h + 
i-o 

(3.19) T 5 io N [q T A v q' - q 1 A„q'_\ + T 7 uj M [q-A v q'_ - qlA v q'] + 

TQio N [q l A„q - q 1 A v q„] + t s lo m [qLA„q- - qiA v q], 


where q 1 denotes the derivative of q either in x or y direction, to is the Legendre weight, and A v is 


(3.20) 


A v - 


/ dFf dG v \ _ 
V 0 q x ’ dq u ) 


Since all the eigenvalues of A,, are non-negative, every term inside the summations in the above equation is 
not negative, and we would like to keep the boundary terms. Thus we get the energy estimate such as 


(3.21) 


< W A v q' -to N q' T A u q' - q f _)A v q'_ - io M q'I A v q'_] + 
T 5 co N [q T A v q' - q T A v q’_] + T 7 io M \q T ~ A v q[__ - q 7 A v q '] + 
T 6 CO N [q T A v q - q T A v q_] + T S ( 0 M [q-A„q~ - q 7 A v q] . 


The RHS of (3.21) can be rewritten as 


(3.22) 


RHS = u 1 Bu, 


where u = (q,q~,q' ,qL) and B is given by 


/ 


(3.23) 


B 


2<T6 A v 
- oyA 7 - o 8 A v 
(1 + 

-cr 5 AT 


— Of, Ay — o 8 A v 

2 0'S Ay 
—a? A 7 
( — 1 + <77) A 7 


(1 + 05 ) A v -or, A v \ 

-07 A v (-1 + ar 7 )A v 

— 2bJpfA v 0 

0 —2 ujmA v J 
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with 0 = diag( 0, 0, 0, 0), <r 5 = u n ts,<j$ = wjv 7e, <r 7 = wa/7y and <rg = wmTs- It is sufficient for the proof if B 
can be shown to be negative semi-definite. This first leads to: 


( 3 . 24 ) wjyre < 0 . uJn'iq = 1 + ukts — — 0 . 

Note that we use here the fact that A„ is symmetrizable (see Appendix B). Taking into account (3.24), B 
becomes 

/ 2 a 6 -a 7 - 1 + er 7 \ 

(3.25) B = — cr 7 -2 wjv 0 

\ — 1 T < r 7 0 — -com ) 


To ensure negative semi-definiteness, det(-B) < 0 and therefore 


1 J <r| — 2 <T 7 + 4(7g 4 — - — < 0. 

% CU/V / WjVf O/'A/ 


(3.27) er 7 < a 7 < erf 

where 

± __ Q-'n ± / (^A r ^’.vf)(I — -'.t/ + iajy)7 

7 WjVf T CJjv y (u.' J Y+W J Vf) 2 

Here we note that the condition that < Tr ■>- must be also satisfied in order for a 7 to have real 
root. This yields the conditions in the Theorem. □ Note that these conditions are given independently of the 
local flow properties. And moreover, the penalty parameters of each domain are constrained by its adjacent 
domain. 

Remark 2. For ft to be outward normal vector the condition (3.17) is now given by 
CO NT® S fl, CON Tq + com t~8 — 0, 1 4" CO A \Tr> 4- co Mir 7 — 0, 


( f- — — ] (J 2 M T$ + 2t 7 -r 4cu ,y Tq 4* T 0 

\COM CON J COM 


with the conditions (3.13) 


3.2.4. The averaging method. We show in this Section that the averaging method can also be written 
as a penalty method with a particular choice of the parameters. 

Euler Equations :. We start first with the Euler equations: consider the following penalty method: 

dq dF dG 

dt + ~dx + ~dy = Tl - 3 G(*>tf)[/ ^ ~ f ( q -)1 + 


where 


T2,iQ( x ,y)[f ( 9 )-/' («-)], 


/'* = (• A ± q x ,A ± q y ) • n\ go , 


Note that the penalty terms use the derivative of the fluxes. 
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Theorem 3.3. // n = 73 — 73 = 74 = | then the above penalty method (3.29) is equivalent to the 

averaging method and is stable. 

Proof. We prove the theorem at the interface x — 0 with the rectangular domain and assume that 

74 


N ----- M . If n = r 3 = 4, and t-> 
(3.31) 


| then the method becomes 


dq 1 1 
St 1 


a? 


n 


dt 


dF 1 d F n 
~8x + dx 


ac 

ay 


and this is obviously equivalent to the Averaging Method. Here note that Following the 

same procedure in Theorem 3.2, the energy equation becomes 

2^7t m = _ 2wv {qlAqI ~ qIlAqH) ,0 + (nqI ~ T:iqII)A+{q * ~ 

(3.32) +(r 2 q 1 - T 4 q II )A~(q I x - q 1 /) | 0 . 

Since n = r 3 — r 2 = r,j = |, and q r (Q.y,t) = q n (Q,y,t), the RHS of the above equation vanishes and 
the energy is bounded by the initial energy. □ 

The Navier- Stokes Equations The averaging method for the N-S equations can be presented as 


(3.33) 


dq dF dG dF,, , dG v 

dt + dx + dy dx + dy 

Ti,zQ{x,y)[f' + {q) - / ,+ (g_)] + 
t2aQ{xaj)[I'~ { q) - /'”(?-)] + 

7"r>,7 Q(x, y ) [A,. • a 2 q - A„ • a 2 q_] + 
n 1 .?,Q{x,y)[A v ■ aq - A,, ■ 5q_], 


where a 2 q is the second derivative of q in either x or y direct ion. 

Theorem 3.4. Ifn = r 3 = §, = 74 = 75 = 77 = end rg = --rg = — 5 ™, then the approximation 

is continuous at the interface and the scheme (3.33) is stable. 

Proof If n = 73 = T‘> = 74 = 75 = 77 = ~. and = - 7 $ = then (3.33) becomes, 

V , dq 11 , __i/ar^ aF J/ \ ar; 

at * _0 at ,r ~° 2 v a® + ax / ay 

1 faiy 1 aF,( 7 \ at;, 

2 \ ax ax / ay 

(3.34) +— A„ • (Oil 11 - aq 7 ), 

OJn 

and this ensures the continuity of the approximation at the interface. If the approximation is smooth enough 
such that the derivative of q is continuous at the interface then this becomes the averaging method. 

Thus we get for the energy: 


(3.35) 


1 d 

2u7/y d t 


E(t) = 


2c djs 


W Aq 1 - q !, Aq n - 2q J A v ql + 2q n ArfJ ) U =0 


-0 f‘2 

q x A v q x dx — / q'J A „q ! J dx 
-2 ' ' -Jo 


+ 


[(nq 1 - nq n )A + + (r 2 q 1 - T 4 q rI )A ] {q\. - qi r )\ x =o + 
[inq 1 - T 7 q n )A v {(£ x - q") + (r 6 g J - r 8 g II )A, / (q I x - q r J)] 


x=0 • 
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Since q 1 (0, y. t) = q n (0, y, t), we have 


| m < /([(T i - r 3 )d + + (75 - n)A-]{ql - ql 1 ) + 

(n - T 7 )A v (cfl xx - qii) 

(3-36) (r 6 - r 8 + — )A„(g* - ))U=o- 

WjV 

Thus if n = rg = |, T 2 = T 4 = To = T 7 = |. and rg = — r 8 = the RHS vanishes. □ 

3.2,5. Adaptive averaging. To ensure the stability of the scheme at some particular collocation points 
where the solution become singular and unstable, vve use the averaging method adaptively at selective grid 
points. In particular, we switched from the penalty method to the averaging when the following criteria was 
satisfied: 


(3.37) 
or 

(3.38) 


m.ax 


\P~P-\ \ T-T.\ 

|p + p_|’|r + r_| 


> C a 


\ P ~ P ~\ 

\P+P-\ 


>C a 


where C ave is a non-negative constant. Note that C ave — 0 leads to the averaging method, whereas a large 
C ave results in the penalty method. For the value of C aue used in this paper, we found out that there were 
very few points in which one needs to switch from the penalty to the averaging procedure. Moreover this 
happened only at very few time steps. 


4. The Cavity System And Numerical Configurations. In this section we describe the set up of 
the simulations of the recessed cavity flameholders by the spectral multi-domain technique presented above. 
The main goal of this experiment is to investigate how the geometry of the aft wall affects the flame stability. 

4.1. Physical setup. In the SCRAMJet community, a cavity with the length-to-depth ratio L/D < 
7 ~ 10 is usually categorized as an ’open’ cavity since the upper shear layer re-attaches at the back face [2]. 
In this work, we choose the L/D of the baseline cavity to be 4 and thus the open cavity system is considered. 
The coordinates of the cavity are (7cm, —lcm) for the upper left and (11cm, —2cm) for the right bottom 
corners of cavity. With the length of the neck of the cavity fixed to be 4cm, we consider three different 
angles of the right corner of the floor of the cavity ( 60,45 and 30), we then compare each one with the 
case of the rectangular aft wall. The fluid conditions are given as followings; the free stream Mach number 
M = 1.91, total pressure P = 2.82 (atm), total temperature T — 830.6(A) and normalized Reynolds number 
R e = 3.9 x 10' (1/m). Note that the Reynolds number is here normalized and has a unit of 1 /[length], 
also the Reynolds number based on the cavity dimensions is O(10 5 ). The boundary layer thickness scale is 
S — 5 x 10 - 4 (m), and finally, the wall temperature is T w — 460.7835(A). The initial configuration for the 
baseline cavity system is shown in Figure 4.1. 

4.2, Numerical setup. We have conducted two different experiments for each of the following cases 
(1) non-reacting cold flow and (2) reacting flow . We use 9 and 17 subdomains for both cases 1 and 2. For 
the outflow conditions at the exit of the system arid at the upper boundary, we mainly use a semi-infinite 
mapping in order to reduce the possible reflections at the boundaries. The characteristic boundary conditions 
are also applied and will be discussed in the next section and compared to the mapping. For the case of the 
reactive flows, the cavity was initially filled with Hydrogen fuel with fuel-to-total gas ratio of 0.5. The order 
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M = IS>I 

Re = 3S>5 E+07 (1M) 
Pt = 2.828522 (atm) 
Tt = 830.6 (K) 

UD = 4 
6= 0.0013 



(d=2mm) 


K 


Fig. 4.1. The initial configuration for the baseline cavity system. 

of the polynomial of approximation in y direction in the domain beside the wall is taken large enough to 
resolve the boundary layer well. Finally the adaptive filtering is turned on if the mass fraction of Hydrogen 
arid Oxygen exceed the range of —0.09 < fn 2 < 1.09. —0.02 < fo 2 < 0.25 and the temperature exceeds the 
range of 300 (A') < T < 3500(A). As the shear layer and the complex features of the flows develop, the 
adaptivity criteria for applying the local smoothing is satisfied at some points. In the calculations, we use 
the 3rd and 2nd order local filtering for the non-reactive and reactive flows respectively. It turns out that 
the local smoothing was applied in very few points at the upper corner of the cavity wall. 

For the adaptive averaging, we use the criteria constant C ave such that the difference of the state vectors 
(or pressure) between the two adjacent domains is less than 10%. In figure 4.2 the Penalty Navier-Stokes 
equations were considered for the non-reactive cold flows. As evident from the contours of the density, the 
approximations were well matched at the interfaces. Here the outer boundary was approximated by using the 
characteristic conditions of the inviscid fluxes. The adaptive averaging, with the given adaptivity conditions 
above, took place at only a few points. The characteristic boundary conditions using the inviscid fluxes yield 
good results for both the problems of the density peak propagation and the non-reactive cold flows. As in 
figure 1. we observe that there exist penalty parameters satisfying the stability conditions that, may induce 
reflecting modes at the interfaces. 


1 















, », 



~A~ 



Fig. 4.2. The non-reactive cold flows with the penalty Navier-Stokes equations: the density contours are given in this 
figure at t — 0.25ms. 17 domains are used and the boundaries of each domain are shown. 


5. Results And Discussion. 

5.1. Pressure history. Figure 5.1 show's the pressure history of the non-reactive cold flows for the 
various angles of the aft wall at two different locations inside the cavity, i.e. at the center, (x, y) = 
(8.5 cm, —1.5 cm), and at the middle of the floor (x,y) = (8.5cm, — 1.9cm). 

These figures show' that the pressure fluctuations in cavities with lower angle of the aft are weaker than 
in cavities with higher angles. It is also shown that the attenuation of the pressure fluctuations are obtained 
both at the center and the middle of the floor of the cavity. It is interesting to observe that the patterns of 
the pressure fluctuations for a given angle at different locations are different depending on the angle. In the 
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Fig. 5.1. Pressure history for non-reactive flows: the left panel represents the pressure history at the center of the cavity 
and the right panel at the middle of the floor of the cavity. Each panel shows the case of 90, 60, 40 and 30 degree cavity walls 
from top to bottom . 



Fig. 5.2. Pressure history of the non-reactive flows with the use of the 4th order filter: the left panel represents the 
pressure history at the center of cavity and the right panel shows the left panel in a smaller scale. Each panel shows the case 
of 90, and 30 degree cavity walls from top to bottom. Note that the scale of the right panel is different from the left. 


case of the 30 degree aft wall, the pressure fluctuations are almost the same at the two locations considered 
whereas the case of 45 degree shows a difference in the patterns of the pressure fluctuations between the two 
locations. The pressure fluctuations at the bottom grows greater than that at the center after some time. 



Fig. 5.3. Stream, lines: the left figure shows the streamlines at t — 1.685ms for the global filtering order 7 = 4 and the 
right at t — 2.38ms for 7 — 6. 

Figure 5.2 shows the pressure history when the heavy global filter is applied (iri this case, the 4th order 
filter was used). Unlike the previous case illustrated in Figure 5.1, where the 6th order global filter is used, 
the pressure fluctuations eventually decay out and a large recirculation zone is formed inside the cavity 
without any severe pressure fluctuations. Note that the scale in the left panel shown is the same as in Figure 
5.1 while the right panel is shown in a smaller scale for a closer look. This figure shows that the large 
recirculation zone(s) formed inside the cavity obtained by the lower order numerical scheme is induced not 
physically but rather artificially due to the heavy numerical dissipations. This is clearly shown in Figure 
5.3. In this figure a large recirculation zone is observed - this zone is formed earlier than this streamlines 
are captured - when the 4th order filter is used (left figure) and an almost steady state is already reached as 
the pressure history indicates in Figure 5.2. We find from the numerical results that the large recirculation 
is very stable once it forms. This large recirculation and the steady state solutions are not observed in the 
case of 7 = 6(right). For the case of 7 = 6 instead of the large single recirculation zone, smaller scale vortex 
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circulations are formed and they are interacting with each other, never reaching the steady state with time. 
This result shows that for these sensitive problems, high order accuracy should be used in order to minimize 
the effect of the numerical dissipation. 

Figure 5.4 shows the case of the reactive flows for the 90 and 30 degree aft walls. Similar features of 
the pressure fluctuations are shown as in the non-reactive flow's. However the pressure fluctuations are much 
more attenuated for both the 90 and 30 degree walls than in the non-reactive cold flows. In the reactive 
cases Hydrogen fuel, which was initially supplied inside the cavity was consumed. As time elapses, the fuel 
is consumed out with the production of the water for these cases. 






Fig. 5.4. Pressure history for reactive flows: the left panel represents the pressure history at the center of cavity and the 
right panel at the middle of the floor of cavity. Each panel shows the case of 90 and 30 degree cavity walls from , top to bottom. 



Fig. 5.5. 


The water contour of the reactive flows: the left the water density contour is given in the left figure and its 


streamlines in the right figure at t = 0.135ms. 


These results demonstrate that simulations of cold flows do not necessarily shed light on the behavior 
of reactive flows. 

5.2. Flow fields. 

Non-reactive cold flow 

Figure 5.6 shows the density contours and streamlines for the 90, 60, 45 and 30 degree walls at the 
instant time t = 2.4ms. As shown in the figure, the shear layer is becoming weaker as the degree of angle 
of the aft wall and the flow fields are becoming more regularized for the case of the lower angle. And note 
that the density compression at the corner of the aft wall is also becoming weaker for the more slanted wall 
cases. 

Figure 5.8 shows the streamlines corresponding to the each case of Figure 5.7. Note that compared to 
the non-reactive cases, the shear layers are less developed for the reactive cases. As the figures of the pressure 
fluctuation history and Figure 5.8 indicate, the shear layers are weak for both the 90 and the 30 degree walls 
in the reactive cases. 

Reactive flow 

Figure 5.7 shows the water contour inside the cavity for the different angles at different time. Here we 
define the region where the flames are generated to be same as the region where the water is produced. As 
the Hydrogen fuel is consumed, the water is produced and starts to be expelled from the cavity to the main 
channel. The flame-holding efficiency is enhanced if the chemical radicals (water in this case) are stably 
circulating and long lasting before they are expelled from the cavity. Figure 5.7 shows that the lower angled 
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Fig. 5.6. The density contour and the streamline of the non-reactive flows: the left column shows the density contour for 
90, 60, f5 and 30 degree walls from top to bottom and the right column shows the corresponding streamlines at t = 2.43 ms. 
The maximum contour level is 1.8 and the. minimum 0.5 with the level step size 50. 


aft wall (30 degree in this case) maintains more water than the 90 degree wall at a given time. The figure 
also shows that the lower angled aft wall holds the flame (water in this case) longer than the 90 degree wall 
- in the last figure in Figure 5.7 at t = 2.26ms, the most water is expelled and the only the small amount 
is left in the left corner while the 30 degree wall cavity holds the water still throughout the cavity. These 
results imply that the flame-holding efficiency can be increased by lowering the angle of the aft wall of the 
cavity. 

Appendix A. The similarity transform matrices and the eigenvalues of the inviscid flux 
with chemical species. 

Air model without combustion 

First consider the ideal gas composed of two chemically non-reactive species (for the ideal mono-atomic 
gas A the diagonal matrix and S, the diagonalizer, were given in [15]). A is given by 

A = diag(U ■ n + c, U ■ ft, U - ft, U - ft — c, U - ft, U ■ ft), 

where U — (u, v ), ft = (n x ,n y ) is an unit outward normal vector at the interface and c is a local sound speed. 
For simplicity we assume that 

p j Cp(s)ds — I’ ~ pC„T , 

Jo 


1(5 




^ = 0.175ms 




t — 0.945ms 



t ------ 2.26ms 


Fig. 5.7. The water contour of the reactive flows: the water density contours are given in the left figures for 90 degree wall 
and 30 degree wail in the right figures. From top to bottom the instant times t are 0.175 ms, 0.275 ms, 0.945 ms and 2.26ms. 
The maximum, and minim, um contour levels are 0.01 and 0.23 respectively with the number of levels 50. 


This form is used only in the analysis, as mentioned iri Sec. 3.1, C Vi is expressed as a 4th order polynomial 
in the temperature T. The nonlinear expression of C Pj makes it difficult to derive the Jacobian matrices 
of the fluxes. Our simplifications is a results of assuming small coefficients of the high order terms of the 
polynomial. Iri the actual simulations C„ is computed appropriately using the empirical law and assumed 
temperature independent at each linearization step. With this assumption S is given by 
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^ = 0.175ms 



t = 0.275ms 



t = 0.945ms 



t — 2.26ms 


Fig. 5.8. The streamlines for the reactive flows: the streamlines for 90 degree wall are shown in the left figures and the 
30 degree wall in the right. From top to bottom the times t are 0.175 ms, 0.275 ms, 0.945ms and 2.51 ms. 


tangential vector k = (— %, n x ) and 


h°; 


(li} ~ Rr (n,h" ' Hjh'iy " ji n ‘ 


h? 


Note that R v = 7 - 1 for the mono-atomic ideal gas with 7 , the ratio between the heat capacities C v and 


Air model with combustion 

Consider now the equations the Euler equations with four reactive species. In this case A and S are 
given by 

A = diag{U ■ n + c, U - ft, U ■ n, V ■ n — c, U ■ ft, U - ft, U ■ ft, U ■ ft), 
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where all the variables are same as in the two species case except that 

dijki = eijki (hj - h\ + h^)R v /R h , 

and 

Rijkl = ~~ ^ijkl ( Rj Rk " t ” Rl)cf Rfii 

with 

4 

Rh - y e m Rdtf - hi + hi), i,j, k.l - 1 . 2, 3, 4, j < k < l, 

i - 1 

eijki is the permutation symbol and R„ = u fa Ri- A and S are based on the time dependent local spatial 
quantities at a given time. f £ is calculated at the interface points at each time. 

Appendix B. The symmetrizability of the coefficient matrices of the Navier-Stokes equa- 
tions with chemical species. 

In [1] it had been proven that the coefficient matrices of the Navier-Stokes equations (expressed in 
the primitive form), of the ideal gas can be simultaneously symmetrized. I 11 [12, 15] the same result was 
demonstrated for the conservative form of the equations. Here we show that it is also true for the Navier- 
Stokes equations of the combustible gas with multiple chemical species in two dimension. 

Rewrite the linearized Navier-Stokes equations (2.2) in conservative form without the chemical source 
term as 

9Q , dq Oq = 0 2 q d 2 q & 2 q 

dt dx dy ' dx 2 dxdy J dy 2 5 

where A=^-,B = ^.C = and E = Of*-. It is sufficient to consider the chemically 

interacting two chemical species. The coefficient matrices are given by 

/ 0 1 0 0 0 0 \ 

X - u 2 (2 - R v )u -R v v R v ip i ip 2 

—uv v u 0 0 0 

A = 

u(x — H) H — R v u~ — RyUV (1 + R„)u uipi uip 2 
-ufi fi 0 0 u 0 

^ -M/2 h 0 0 0 u j 
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where U 2 = \U U,H = 
R v h°i,Si = -hifi^Ji = 


E+P 
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o 3 ’ 
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= ^,02 = %,o* = M^r,C = jfc.X = «,(C'2 - C),% = CRi 
f + E;=1 ( ft * - $i)fi,0 = £~ 2 f / 2 and £>•* = <t,u 2 + cr k v 2 + a 3 0 . 


To find the symmetrizer for A, B.C.D and E, we first consider the similarity transform matrix S F of 
C such that 


S^CS p = Ac, 


where Ac is a diagonal matrix composed of the eigenvalues of C. The subscript P denotes that this matrix 


adopted from the parabolic portion of the equations [1] 

. The diagonal matrix Ac of C is given by 
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Introducing a symmetrizing diagonal matrix, Q T Q such as 



The diagonalizer S P is composed of the eigenvectors of C, S P and its inverse S p 2 are given by 



we have symmetrized all the coefficient matrices, i.e. 


Q T QSp 1 AS P = (Q T QSp 1 AS P ) T , 

Q T QSp 1 BS P = (Q T QSp 1 BSp) T , 

Q T QSp 1 CSp - (Q T QSp 1 CSp) T , 

Q T QSp 1 DSp - (Q T QSp 1 DSp) T , 

Q T QSp 1 ESp = Sp 1 ES P . 

Appendix C. Constants for Chemical Models. 

Here we provide constants used in the chemical model for the current numerical experiment. Table I gives 
the constants used to get the approximation of the specific heat C p i of i th species in the 4th order polynomial 
of T. i.e. 


Cpi = (ci + T(c 2 + T(c 3 4- T(c 4 + csT))))R/M i} 
where R is a gas constant, and Mj is a molecular weight of i th species [16]. 


Table I 

Coefficients for the approximation of the specific heat C p i 



0 2 

Hi 

H 2 0 

N 2 

ci ( 1 / mole) 

3.0809 

3.4990 

3.4990 

3.1459 

02 ( 1 /mole) 

0.16962E-2 

-0.18651E-3 

0.14878E-2 

0.99154E-3 

cs(l /mole) 

-0.76334E-6 

0.46064E-6 

0. 87544 E-7 

-0.22912E-6 

04(1 /mole) 

0.17140E-9 

-0.13157E-9 

-0.11499E-9 

0.12181E-10 

05 ( 1 /mole) 

-0.14116E-13 

0.11679E-13 

0.1.3495E-13 

0.11024E-14 


Table II gives the molecular weight and specific enthalpy for each chemical species and Table III gives 
the reference dynamic viscosity, temperature constants T and S in Wilke’s law [21]. 


Table II 

Molecular weights and specific enthalpy 




0-2 

Hi 

h 2 o 

N 2 

M (1 /mole) 


32.000 

2.016 

18.016 

28.016 

hP (Joule/ kg) 

-272918.21 

-4280070.46 

-13973684.55 

-302736.23 

Table III 

Constants for Wilke’s law 


0-2 

Hi 

H 2 0 

iV 2 

fio(kg jrn / sec) 

0.1919E-4 

0.0841 IE- 4 

0.1703E-4 

0.1663E-4 

T 0 (K) 

273.111 

273.111 

416.667 

273.111 

S(K) 

138.889 

96.6667 

861.111 

106.667 
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